% 插值多项式
% Gitee Repo

clc
clear

global N1
global N2
global u

% f(x) = u(1) + u(2)x + u(3)x^2+...+u(N1+N2)x^(N1+N2-1)
% (N1+N2) 个点插值 (N1+N2-1)阶多项式
% N1 (x1,y1): 函数值，即 y1=f(x1)
% N2 (x2,y2): 函数导数值，即 y2=f’(x2)

N1 = 3;
%x1,y1:(N1,1)
x1 = [0;1/2;1];
y1 = [0;1;-1/2];

N2 = 3;
%x2,y2:(N2,1)
x2 = [0;1/2;1];
y2 = [0;0;0];

% check
if ndims(x1)!=2 || ndims(y1)!=2 ...
    || size(x1,1) != N1 || size(y1,1) != N1 ...
    || size(x1,2) != 1 || size(y1,2) != 1
    error('x1 dim mismatched')
end
if N2 > 0
    if ndims(x2)!=2 || ndims(y2)!=2 ...
        || size(x2,1) != N2 || size(y2,1) != N2 ...
        || size(x2,2) != 1 || size(y2,2) != 1
        error('x2 dim mismatched')
    end
end

A = zeros(N1+N2,N1+N2);
b = zeros(N1+N2,1);

A1 = x1.^(0:(N1+N2-1));
A1(:,1) = 1;

if N2 > 0
    A2 = (0:(N1+N2-1)).*x2.^((0:(N1+N2-1))-1);
    A2(:,1) = 0;
    A2(:,2) = 1;
end

if N2 > 0
    A = [A1;A2];
    b = [y1;y2];
else
    A = A1;
    b = y1;
end

% u: (N1+N2,1)
u = A\b;

function res = f(x)
    global N1
    global N2
    global u

    if ndims(x) != 2 || size(x,2) != 1
        error('x should be (N,1)')
    end

    res = x.^(0:(N1+N2-1));
    res(:,1) = 1;
    res = sum(u'.*res,2);

    return;
end

figure
hold on
#axis equal

x = (0:0.01:1)';
plot(x,f(x));

% 绘制点
scatter(x1,y1)

% 绘制导数点
if N2 > 0
    f2 = f(x2);
    scatter(x2,f2)
    for i = 1:N2
        line([x2(i)-0.05,x2(i)+0.05],[f2(i)-0.05*y2(i),f2(i)+0.05*y2(i)])
    end
end

